function u = df_solve(N,P,Tree,f)
  
  L = round(log(N/P)/log(2))+1;

  u = f;
  %-----------------
  u = f;
  for ell=1:L
    W = 2^(ell-1)*P;
    nck = N/W;
    for i=1:nck
      for j=1:nck
        TT = Tree{ell}{i,j};        %LINS = TT.LINS;        LBDS = TT.LBDS;
        GINS = TT.GINS;        GBDS = TT.GBDS;
        invA = TT.invA;        B = TT.B;        S = TT.S;
        %load
        in = u(GINS,:);
        bd = u(GBDS,:);
        tt = invA*in;
        %L
        bd = bd - B*tt;
        %solve
        in = tt;
        if(ell==L)
          bd = S*bd;
        end
        %save
        u(GINS,:) = in;
        u(GBDS,:) = bd;
      end
    end
  end
  %-----------------
  for ell=L:-1:1
    W = 2^(ell-1)*P;
    nck = N/W;
    for i=1:nck
      for j=1:nck
        TT = Tree{ell}{i,j};        %LINS = TT.LINS;        LBDS = TT.LBDS;
        GINS = TT.GINS;        GBDS = TT.GBDS;
        invA = TT.invA;        B = TT.B;
        %load
        in = u(GINS,:);
        bd = u(GBDS,:);
        %Lt
        in = in - invA*(B.'*bd);
        %save
        u(GINS,:) = in;
        u(GBDS,:) = bd;
      end
    end
  end
  %-----------------
  
  
